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Manned Spacecraft Center 


SUMMARY 


The primary concern of this paper is to investigate the problem of inver- 
sion of singular or nonsquare matrices . In this connection, a new algorithm 
for computing the generalized inverse of an arbitrary complex matrix is given. 
For a nonsingular matrix the algorithm gives the ordinary inverse of the matrix ■ 

The paper is divided into several sections. The first two sections give 
a definition- theorem expose' of the known results in the literature. The fol- 
lowing sections give a new explicit form, together with an algorithm for com- 
puting the new explicit form. An application to least squares approximation 
is given that can easily be realized in trajectory analysis problems. Finally, 
a computer program for computing the generalized inverse of a matrix is given 
utilizing the algorithm mentioned in the latter paragraph. 


INTRODUCTION 


A. Bjerhammar (ref. l), E. H. Moore (ref. 2), and R. Penrose (ref. j) 
independently generalized the concept of matrix inversion to include arbitrary 
complex matrices. Their equivalent forms of the generalized inverse of a 
matrix have given rise to many applications of generalized inversion. In the 
text, the basic theory is utilized in giving a new explicit form of the gener- 
alized inverse of an arbitrary complex matrix. 


SYMBOLS 


Capital letters 
Lowercase letters 

A* 

A ” 1 


matrices unless otherwise stated 

column vectors unless otherwise stated or clear from 
context 

matrix conjugate transpose of A 
matrix inverse for nonsingular A 



generalized inverse of A 


A 

H 

R(A) 

P R(A) 

E 111 


diag^a 1 ,a 2 , . . .,a^ 


a Hermitian idempotent matrix (h.i.); that is, a matrix 
such that E* = H and HE = H 

range space of A ; that is, the collection of all images 
of column vectors under the transformation A 

orthogonal projection on the range of A; that is, a 
Hermitian idempotent leaving R(A) fixed 

m- dimensional euclidean space 

diagonal matrix 


DEFINITIONS AND EQUIVALENT FORMS 


A. B j erhammar (ref. l), E. H. Moore (ref. 2), and B. Penrose (ref. 3) 
independently generalized the concept of matrix inversion^ to include arbitrary 
complex matrices. The generalized inverse of a singular, or nonsquare, matrix 
possesses properties which make it a central concept in matrix theory. 

In this, paper, a definition- theorem expose is presented, along with 
applicable references and special problems. The following fundamental theorem 
due to Penrose (ref. 3) is stated without proof: 

THEOREM I. The four equations 


AXA = A 

(1) 

XAX = X 

(2) 

(AX)* = AX 

(3) 

(XA)* = XA 

W 


have a unique solution X for each complex matrix A. 

Definition 1. The solution X in THEOREM 1 is denoted X = A + and is 
called the generalized inverse of A. 

The following theorem gives an equivalent form of A . 

THEOREM II . For any mxm matrix A over the complex field, X = A + 
is the unique solution to the equations AX = P R ( A ) and ^ = P r(x) 
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is the orthogonal 


where R(A) is the range space of A in E 131 and P R(A) 
projection on R(A) . 

Proof: THEOREM I implies that AX is a Hermitian idempotent (see Symbols) 

which leaves A fixed, that is, (AX)A = A. Hence, AX must be a projection. 
It may be concluded that XA is also a projection. 

The properties of the generalized inverse and possible computing schemes 
are. given in the following theorem. 

THEOREM xil . Let A be an arbitrary complex matrix. Then, for scalar 
\ / 0 and unitary U and V, 


+ . +.* » + +.* + 
A + (A ) A = A = A (A ) A 

(5) 

+ # * # + 
A AA = A = A AA 

(6) 

(A + ) + = A 

(T) 

<aV > (aY 

(8) 

A + = A" 1 for nonsingular A 

(9) 

/ — . 
II 

(10) 

(A*A) + = A + (A + )* 

(11) 

(UAV) + = V _1 A + U _1 

V * ^ 

A = ) A^ and A^Aj = 0 

(12) 

A% t = 0 for i ^ j 

imply A + = A? 

J 

> (13) 
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If A 


is normal (i.e., A A = AA ) then 

A + A = AA + and (A n ) + = (A + ) n 
* + + . 

A, A A, A , and A A all have rank equal to 
trace (A + A) 

+ , * .+ * 

A = (A A) A 


(l4) 


(15) 

( 16 ) 


Equation (l6) reduces the problem of computing A + to that of computing 

the generalized inverse of a Hermitian matrix A A. Moreover, such a matrix 
can always be diagonalized by a unitary transformation, that is, 

D = U(A* A) V = diag^, ...,a n j 

Equations (10) and (12) imply that 


(A*A) + = VD + U = V diag^, . . .,^U 

It is tacitly assumed that if a. = 0, the corresponding term in 
/l 1 \ 1 

diag — , . ..,— I is zero. It is not usually an easy task to determine the 
\ a i a n / 

unitary transformations U and V. Methods for computing the generalized 
inverse have been given by various authors (refs. 1, k, 5, 6, and 7)* 

The following is a theorem of major importance characterizing all sole 
tions of the matrix equations AXB = C which have some solution X. 

THEOREM IV . For the matrix equation AXB = C to have a solution, a 
necessary and sufficient condition is 

AA + CB + B = C 

in which case, the general solution is 

X = A + CB + + Y - A + AYBB + 

where Y is arbitrary to within the limits of being consistent with the di- 
mension in the indicated multiplications (ref. 3) . 
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Proof: If X satisfies AXB = C, 


C = AXB = AA + AXBB + B = AA + CB + B 

_l_ _L + + 

Conversely, if C = AA CB B, A CB is a particular solution. Clearly, for 
the general solution, AXB = 0 mist he solved. Any expression of the form 

X = Y - A + AYEB + 

-ft + 

The only property required of A and B in the theorem is AA A = A and 
EB + B = B. 


Corollary A. The general solution to the vector equation 


Ex = c 


is 


x = P + c + (I - P + P)y 

where y is arbitrary, provided a solution exists. 

Corollary B. A necessary and sufficient condition for the equations 

AX = C 

and 

XB = D 


to have a common solution is that each have a solution and AD = CB (ref. 8). 

Proof: If AX = C and XB = D have a common solution, then clearly each has 

a solution and 

AXB > CB 
AXB = AD 

so that 


CB = AD 

In order to obtain the sufficiency of the condition, it is assumed that 

X = A + C + DB + - A + ADB + 

which is a solution if AD = CB, AA + C = C, and DB + B = D. 
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THEOREM V. 


A + A, AA + , I - A + A, and I - AA + are h.i. (17) 

H is h.i. implies H + = H (l8) 


(See Symbols.) 

Proof: The proof requires a straightforward application of THEOREM I. 

In general, the reversal rule (i.e., (AB) + = B + A + as in the case of the 
standard inverse) does not hold. R. Cline (ref. 9) recently obtained the 
following result: 


THEOREM VI . 

Then, 
where : 


and 


Let A and B be matrices with the produce AB defined. 
(AB) + = B+A* 

AB = AjB x 
B ± = A + AB 


THE EXPLICIT FORM 


Utilizing the properties of A + in the preceding sections, an explicit 
form is developed which gives rise to an algorithm for computing the generalized 
inverse of an arbitrary complex matrix (ref. 5) • 

THEOREM VII . For any matrix A, A + = WAY where W and Y are any 
solutions of 

WAA* = A* (19) 

and 

* * 

A AY = A 

Proof: Equations (19) and (20) indeed have a solution W 

if W and Y are any solutions, 

* * * * 

AWAA = AA and A AYA = A A 


( 20 ) 


= Y = A . Moreover, 


so that 


AWA = A and AYA = A 
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Note: BAA* = CAA* implies BA = CA 


In addition. 


* # * * # * * * 

WAA W = A W and Y A AY = Y A 


imply 

(WA)* = WA and (AY)* = AY 

If it is assumed that X = WAY, X satisfies equations (l) to (4) of 
THEOREM I so that A + = X = WAY. 

+ * * 

Corollary C. For any matrix A, A = A S^AS^A where and are, 

respectively, any solutions of 

(AA r )S 1 (AA ) = (AA* ) 

and 

(A # A)S 2 (A*A) = (A*A) 

Proof: According to THEOREM III, W = A and Y = S^A are solutions of 

equations (19) and (20) of THEOREM VII provided 


(AA*)S l (AA*) = (AA*) 

and 

(A A)S 2 (A A) = (A A) 

The corollary follows . 

THEOREM VTII , If B is a matrix and nonsingular matrices P and Q 
exist so that FBQ = E is an idempotent, then B = QEP is a solution of 
BXB = B. 

Proof: If P, Q, and E satisfy the hypothesis of the theorem. 


B = P" 1 EQ" 1 

and 


BBB = ^p" 1 EQ,“ 1 ^ QEP^P" 1 EQ” 1 j = P -1 EQ -1 = B 

Corollary C and THEOREM VIII suggest an algorithm for computing the 
generalized inverse of a complex matrix F. Consider the equation 
F* =* (F F) + F from reference 10, which reduces the problem of finding F 
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* 

to that of finding the generalized inverse of the Hermitian matrix F F = C. 

Since (c 2 )* = C 2 , nonsingular matrices P and Q exist (products of elemen 
tary matrices obtained by simple elimination) so that 


pc 2 q = 



where I is a rank r identity matrix and the Z is a zero matrix. When 
C is set equal to A in Corollary C, 


* * * # 
A A = AA = C C = CC 


2 


C 


According to THEOREM VIII, solutions S_, = S 0 = QI P are chosen so that 

1 £ o 


C + = (CS 1 ) 2 C 


and finally. 


(F*F) + = C + 
F + = C + F* 


Computing programs for calculating and are now in existence 

(e.g., STORM, Statistically Oriented Matrix Program, IBM). In general, these 
programs only compute some solution of the equation AXA = A, usually different 
from A + . These results allow one to construct a solution to all four Penrose 
equations (eqs. (l) to (4) of THEOREM I), given only a solution of the first, 
namely, AXA = A. 


APPLICATION TO LEAST SQUARES APPROXIMATION 


In this section, an application to the least squares approximation is 
stated that can be realized in trajectory analysis problems. For the sake of 
simplicity, weighting is not considered; however, it would introduce no 
difficulty. 

The vector equation Ax = b does not, in general, have a solution x. 
However, all candidates for a least squares solution (i.e., a solution vector 
x minimizing (Ax - b) (Ax - b)) must be solutions of the normal equations 
(ref. 6). 


* * 

A Ax = A b 
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THEOREM IX . 
The equation 


Let A 


be any matrix 


and b 


be any vector (mxl) . 


(mxn) 


A*Ax = A*b 


always has a solution, and hence a general solution is given by 

X = (A A.) + A b + [i - (A*A) + AVjy 
■ A + b + (I - A + A)y 

* 

Moreover, if A A is nonsingular, the solution is 

x = A + b 


and is unique. 

Proof: First, it is shown that 

A Ax = A b (21) 

has a solution. Consider the vector 

x = A + b 

* + v * + 

Since equation (6) of THEOREM III implies A A(A b) - A b, x = A b is indeed 
a solution of equation (21). The existence of this solution, together with 
Corollary A, implies that the general solution to equation (21) is 

x = (A*A) + A*b + [i - (A*A) + A% *jy (22) 

By using equation (l6) of THEOREM III, 

x = A + b + (I - A + A)y 

* 

Finally, if A A is nonsingular, then 

x = (A A; A b + (I - l)y 
= A + b 

and equation (21) has a unique solution. 

In summary, if x is a least squares solution of Ax = b, then x must 
satisfy 

«• # 

A Ax - A b 
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All solutions of this equation are given by 

x = A + b + (i - A + A)y 

Any vector of the form 

x = A + b + (I - A + A)y 

is a "candidate" for a least squares solution, and this form describes the 
"class of all candidates." 

Corollary D. Every solution of 

-if- 

A Ax = A b 

minimizes 

Q = (Ax - b) (Ax - b) 
provided Q has a minimum. 

Proof: Any vector at which Q is minimum is of the form 

x = A + b + (I - A + A)y 

If Q has a minimum, let 

x 1 = A + b + (I - A + A)y 2 

be any other solution. In order to show that 

(Ax-j^ - b) (Ax 1 - b) = (Ax 2 - b)* (Ax g - b) 

Ax^ and Ax 2 is examined in light of THEOREM I. 

AXp = A^A + b + (I - A + A)yJ 

= AA + b + (A - AA + A)y 1 

= AA + b + (A - A)y 

= AA + b 
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Similarly, 


Axg = AA + b 

so that 

Ax 1 - b)^(Ax 1 - b) = (Ax 2 - b) (Ax 2 - b) 

that is, every vector of the form 

x = A + b + (I - A + A)y 
yields the same minimum value of Q. 

SUBROUTINE GENINV 


GENINV is a FORTRAN IV subroutine used to compute the generalized inverse 
of an mxn matrix A. All computations are made in double- precision floating 
point arithmetic. The subroutine Is based on the algorithm suggested by the 
explicit form. 


Calling Sequence 
Call GENINV (A, AP, M, N, L, E) 
where: 

A double-dimensioned, double- precision array containing the original matrix. 
A is dimensioned A(25, 25 ) 

AP double-dimensioned, double-precision array where the generalized inverse 
of A will be computed. AP is dimensioned AP(25, 25) 

M number of rows in the original matrix 

N number of columns in the original matrix 

L twice N 

E some small number for near-zero divisor test 
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Method 


Given A 
Compute : 


Find nonsingular matrices 


* 

C = A A 

_2 

C = CC 

E and P such that 


(Print A) 

(Print C) 
(Print C 2 ) 



(Print E, P, I q ) 

(A form of Gaussian elimination with pivoting employed) 

Compute: 

R = PI E (Print R) 

o 

then 


C + = CRCRC 


also 


A + = 


+ * 
C A 


(Print C + ) 


(Print A + ) 


Remarks 


The program uses two double- precision arrays CSQ(50, 50 ) and B(25, 25) 
for internal manipulation. The subroutine leaves the original matrix A intact. 

Results are printed after each step as indicated. 


CONCLUSION 


The theory developed in THEOREMS VII and VIII gives rise to easy calcu- 
lation of the generalized inverse of an arbitrary complex matrix when only 
the solution to the matrix equation AXA = A can be found. In general, a 
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simultaneous solution must be found for four matrix equations, 
THEOREM I, that define the generalized inverse. 


given in 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, February 5 , 19 65 
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